note: Using 1-d full PIC code written in FORTRAN (other files).
%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
import math as mt
from scipy.fftpack import fft, ifft,fft2,ifft2,fftshift,ifftshift
from IPython import display
import matplotlib.animation as animation
plt.rcParams['animation.html'] = 'html5'
bs=1024
nb=0.2
#fig,(ax1,ax2)=plt.subplots(1,3,figsize=(10,5))
fig=plt.figure(figsize=(12,4))
ax1=plt.subplot2grid((1,3),(0,0),colspan=2)
ax2=plt.subplot2grid((1,3),(0,2))
def update_anim(i):
it=i*16
f1=open('xp1_%05d' %(it),'rb')
f2=open('xp2_%05d' %(it),'rb')
f3=open('xp3_%05d' %(it),'rb')
xp1=np.fromfile(f1,dtype='float64').reshape(-1,4)
xp2=np.fromfile(f2,dtype='float64').reshape(-1,4)
xp3=np.fromfile(f3,dtype='float64').reshape(-1,4)
for ax in (ax1,ax2):
ax.clear()
ax1.plot(xp1[:,0],xp1[:,1],'C0,')
ax1.plot(xp2[:,0],xp2[:,1],'C1,')
ax1.plot(xp3[:,0],xp3[:,1],'C2,')
ax1.set_xlabel(r'$x/(c/\omega_{pe})$')
ax1.set_ylabel(r'$v_x/c$')
ax1.set_ylim(-0.5,0.5)
ax2.hist(xp1[:,1],bins=bs,range=(-0.25,0.25),histtype='step',color='b',weights=(1-nb)*np.ones(len(xp1)))
ax2.hist(xp2[:,1],bins=bs,range=(-0.25,0.25),histtype='step',color='orange')
ax2.hist(xp3[:,1],bins=bs,range=(-0.25,0.25),histtype='step',color='g',weights=nb*np.ones(len(xp1)))
ax2.set_ylabel(r'$v_x/c$')
ax2.set_xlabel(r'$f(v_x)$')
ax2.set_ylim(0,6000)
anim=animation.FuncAnimation(fig,update_anim,frames=64)
plt.close()
anim
%%time
dx=0.1
dt=0.1
isav=32
nt=4096
nx=1024
wce=0.2
rm=25
for it in range(1,3):
f1=open('ex%03d.hst' %(it),'rb');f2=open('ey%03d.hst' %(it),'rb');f3=open('ez%03d.hst' %(it),'rb')
f4=open('by%03d.hst' %(it),'rb');f5=open('bz%03d.hst' %(it),'rb')
ex=np.fromfile(f1,dtype='float64').reshape(-1,nx)
ey=np.fromfile(f2,dtype='float64').reshape(-1,nx)
ez=np.fromfile(f3,dtype='float64').reshape(-1,nx)
by=np.fromfile(f4,dtype='float64').reshape(-1,nx)
bz=np.fromfile(f5,dtype='float64').reshape(-1,nx)
if it==1:
exsave=ex;eysave=ey;ezsave=ez;bysave=by;bzsave=bz
if it>1:
exsave=np.r_[exsave,ex]
eysave=np.r_[eysave,ey]
ezsave=np.r_[ezsave,ez]
bysave=np.r_[bysave,by]
bzsave=np.r_[bzsave,bz]
%%time
xmax=nx*dx
tmax=nt*isav*dt/rm*wce
fig=plt.figure(figsize=(20,10))
ax1 = fig.add_subplot(231)
ax2 = fig.add_subplot(232, sharex=ax1, sharey=ax1)
ax3 = fig.add_subplot(233, sharex=ax1, sharey=ax1)
ax4 = fig.add_subplot(235, sharex=ax1, sharey=ax1)
ax5 = fig.add_subplot(236, sharex=ax1, sharey=ax1)
im1=ax1.imshow(exsave,aspect='auto',origin='lower',cmap='terrain_r',extent=[0,xmax,0,tmax])
im2=ax2.imshow(eysave,aspect='auto',origin='lower',cmap='bwr',extent=[0,xmax,0,tmax],clim=(-0.01,0.01))
im3=ax3.imshow(ezsave,aspect='auto',origin='lower',cmap='bwr',extent=[0,xmax,0,tmax],clim=(-0.01,0.01))
im4=ax4.imshow(bysave,aspect='auto',origin='lower',cmap='bwr',extent=[0,xmax,0,tmax],clim=(-0.1,0.1))
im5=ax5.imshow(bzsave,aspect='auto',origin='lower',cmap='bwr',extent=[0,xmax,0,tmax],clim=(-0.1,0.1))
ax1.set_xlabel(r'$x/(c/\omega_{e})$');ax4.set_xlabel(r'$x/(c/\omega_{e})$');ax5.set_xlabel(r'$x/(c/\omega_{e})$')
ax1.set_ylabel(r'$t\Omega_{ci}$');ax4.set_ylabel(r'$t\Omega_{ci}$')
ax1.set_title(r'$E_x$');ax2.set_title(r'$E_y$');ax3.set_title(r'$E_z$')
ax4.set_title(r'$B_y$');ax5.set_title(r'$B_z$');
plt.colorbar(im1,ax=ax1)
plt.colorbar(im2,ax=ax2)
plt.colorbar(im3,ax=ax3)
plt.colorbar(im4,ax=ax4)
plt.colorbar(im5,ax=ax5)
plt.show()
nt=4096
enex=np.zeros(nt)
eney=np.zeros(nt)
enez=np.zeros(nt)
enby=np.zeros(nt)
enbz=np.zeros(nt)
for it in range(nt):
enex[it]=np.average(exsave[it,:]**2)
eney[it]=np.average(eysave[it,:]**2)
enez[it]=np.average(ezsave[it,:]**2)
enby[it]=np.average(bysave[it,:]**2)
enbz[it]=np.average(bzsave[it,:]**2)
tt=np.linspace(0,nt*dt/rm*wce,nt)
plt.plot(tt,enex)
plt.plot(tt,eney)
plt.plot(tt,enez)
plt.plot(tt,enby)
plt.plot(tt,enbz)
plt.xlabel(r'$t\Omega_{ci}$')
plt.ylabel(r'$E_x^2$, $E_y^2$, $E_z^2$, $B_y^2$, $B_z^2$')
plt.yscale('log')
plt.ylim(1e-7,1e-2)
plt.show()
%%time
bs=1024
nt=1024
den1save=np.zeros((nt,bs))
den2save=np.zeros((nt,bs))
den3save=np.zeros((nt,bs))
avl1save=np.zeros((nt,bs))
avl2save=np.zeros((nt,bs))
avl3save=np.zeros((nt,bs))
tmp1save=np.zeros((nt,bs))
tmp2save=np.zeros((nt,bs))
tmp3save=np.zeros((nt,bs))
for it in range(nt):
f1=open('xp1_%05d' %(it),'rb')
f2=open('xp2_%05d' %(it),'rb')
f3=open('xp3_%05d' %(it),'rb')
xp1=np.fromfile(f1,dtype='float64').reshape(-1,4)
xp2=np.fromfile(f2,dtype='float64').reshape(-1,4)
xp3=np.fromfile(f3,dtype='float64').reshape(-1,4)
(den1,bns)=np.histogram(xp1[:,0],bins=bs)
(den2,bns)=np.histogram(xp2[:,0],bins=bs)
(den3,bns)=np.histogram(xp3[:,0],bins=bs)
(avl1,bns)=np.histogram(xp1[:,0],bins=bs,weights=xp1[:,1])
(avl2,bns)=np.histogram(xp2[:,0],bins=bs,weights=xp2[:,1])
(avl3,bns)=np.histogram(xp3[:,0],bins=bs,weights=xp3[:,1])
tmp1=0.5*((xp1[:,1]-avl1[np.int_(xp1[:,0])]/den1[np.int_(xp1[:,0])])**2+xp1[:,2]**2+xp1[:,3]**2)
tmp2=0.5*((xp2[:,1]-avl2[np.int_(xp2[:,0])]/den2[np.int_(xp2[:,0])])**2+xp2[:,2]**2+xp2[:,3]**2)
tmp3=0.5*((xp3[:,1]-avl3[np.int_(xp3[:,0])]/den3[np.int_(xp3[:,0])])**2+xp3[:,2]**2+xp3[:,3]**2)
(temp1,bns)=np.histogram(xp1[:,0],bins=bs,weights=tmp1)
(temp2,bns)=np.histogram(xp2[:,0],bins=bs,weights=tmp2)
(temp3,bns)=np.histogram(xp3[:,0],bins=bs,weights=tmp3)
binx=0.5*(bns[1:]+bns[:-1])
den1save[it,:]=den1
den2save[it,:]=den2
den3save[it,:]=den3
avl1save[it,:]=avl1/den1
avl2save[it,:]=avl2/den2
avl3save[it,:]=avl3/den3
tmp1save[it,:]=temp1/den1
tmp2save[it,:]=temp2/den2
tmp3save[it,:]=temp3/den3
%%time
plt.figure(figsize=(20,15))
plt.subplot(331);plt.imshow(den1save,origin='lower',aspect='auto',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.title(r'$density$')
plt.subplot(332);plt.imshow(avl1save,origin='lower',aspect='auto',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.title(r'$mean velocity$')
plt.subplot(333);plt.imshow(tmp1save,aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.title(r'$temperature$')
plt.subplot(334);plt.imshow(den2save,aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.subplot(335);plt.imshow(avl2save,origin='lower',aspect='auto',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.subplot(336);plt.imshow(tmp2save*25,aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.subplot(337);plt.imshow(den3save,aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.xlabel(r'$x/(c/\omega_{e})$')
plt.ylabel(r'$t\Omega_{ci}$')
plt.subplot(338);plt.imshow(avl3save,origin='lower',aspect='auto',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.subplot(339);plt.imshow(tmp3save*25,aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='terrain_r');plt.colorbar()
plt.show()
dx=0.1
dt=dx
isav=32
nx=1024
nt=4096
kmin=2*mt.pi/(dx*nx)*(-nx/2);kmax=2*mt.pi/(dx*nx)*(nx/2)
wmin=2*mt.pi/(dt*nt*isav)*(-nt/2);wmax=2*mt.pi/(dt*nt*isav)*(nt/2)
fftbysave=fftshift(fft2(bysave[:,::-1])) #switching the raw
#plt.pcolor(np.abs(fftexsave),cmap='terrain')
plt.imshow(np.abs(fftbysave),origin='lower',extent=[kmin,kmax,wmin,wmax],aspect='auto',cmap='terrain')
plt.axis([0,kmax/10,0,wmax])
plt.colorbar()
plt.clim(0,100)
plt.xlabel(r'$k(c/\omega_{pe})$')
plt.ylabel(r'$\omega/\omega_{pe}$')
plt.show()
%%time
nt=4096
nx=1024
plt.figure(figsize=(12,8))
fftbz=fftshift(fft2(bzsave[:,:]))
bzpls=np.zeros((nt,nx), dtype=np.complex)
bzmin=np.zeros((nt,nx), dtype=np.complex)
plt.subplot(131)
plt.imshow(bzsave,aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='bwr',clim=(-0.1,0.1))
plt.colorbar()
plt.title(r'$B_z^{raw}$')
plt.xlabel(r'$x/(c/\omega_{e})$')
plt.ylabel(r'$t\Omega_{ci}$')
plt.subplot(132)
bzpls[0:nt//2,nx//2:nx]=fftbz[0:nt//2,nx//2:nx]
bzpls[nt//2:nt,0:nx//2]=fftbz[nt//2:nt,0:nx//2]
bzpls=ifft2(ifftshift(bzpls))
plt.imshow(np.real(bzpls),aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='bwr',clim=(-0.1,0.1))
plt.colorbar()
plt.title(r'$B_z^+$')
plt.subplot(133)
bzmin[0:nt//2,0:nx//2] =fftbz[0:nt//2,0:nx//2]
bzmin[nt//2:nt,nx//2:nx]=fftbz[nt//2:nt,nx//2:nx]
bzmin=ifft2(ifftshift(bzmin))
plt.imshow(np.real(bzmin),aspect='auto',origin='lower',extent=[0,xmax,0,tmax],cmap='bwr',clim=(-0.1,0.1))
plt.colorbar()
plt.title(r'$B_z^-$')
plt.show()
it=250
iphs=128
bs=200
f1=open('xp1_%05d' %(it),'rb')
f2=open('xp2_%05d' %(it),'rb')
f3=open('xp3_%05d' %(it),'rb')
xp1=np.fromfile(f1,dtype='float64').reshape(-1,4)
xp2=np.fromfile(f2,dtype='float64').reshape(-1,4)
xp3=np.fromfile(f3,dtype='float64').reshape(-1,4)
xion=np.r_[xp2,xp3]
plt.hist2d(xion[:,1],xion[:,2],bins=bs,cmap='terrain_r')
plt.xlabel(r'$v_x/c$')
plt.ylabel(r'$v_y/c$')
plt.title(r'$t\Omega_{ci}=%.2f$' %(it*iphs*dt/rm*wce))
plt.show()
!cat input.f90